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ABSTRACT 


A  discussion  is  given  of  certain  methods  of  importance  sampling  and  scoring 
in  the  Monte  Carlo  solution  of  the  radiation  transport  equation. 
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1.  INTRODUCTION 


In  the  spring  of  1961,  the  Development  Division  of  United  Nuclear  Corpora¬ 
tion,  with  the  sponsorship  of  the  U.S.  Army  Nuclear  Defense  Laboratory,  under¬ 
took  a  study  aimed  at  the  preparation  of  a  comprehensive  radiation  transport  pro¬ 
gram.  A  survey  was  made  of  existing  codes  capable  of  yielding  information  on 
the  penetration  of  neutrons  and  gamma  rays  in  the  atmosphere  following  a  nuclear 
weapons  detonation.  In  addition,  information  was  obtained  from  individuals  inter¬ 
ested  in  weapons  effects  on  the  type  and  range  of  answers  needed. 

On  the  basis  of  these  surveys,  a  set  of  recommendations  was  prepared  for 
the  new  program.'  In  addition,  flow  charts  for  portions  of  the  code  were  prepared. 
When,  however,  the  Los  Alamos  Scientific  Laboratory  volunteered  to  prepare  the 
actual  machine  program  for  the  STRETCH  computer,  the  United  Nuclear  activity 
was  redirected  so  as  to  provide  support  to  Los  Alamos.  This  support  has  been 
in  the  areas  of  investigation  of  Monte  Carlo  techniques  and  the  preparation  of 
neutron  cross  sections  for  the  program. 

The  latter  activity  is  reflected  in  the  second  part  of  this  report  (Volume  B) 
which  contains  a  set  of  cross  section  information  for  Be*. 

This  part  of  the  report  is  a  description  of  studies  of  certain  topics  related 
to  the  Monte  Carlo  solution  of  the  transport  equation.  A  really  comprehensive 
solution,  as  now  planned,  will  tax  the  speed  of  the  fastest  present  computer  so 
that  methods  of  improving  the  computational  efficiency  are  needed. 

In  undertaking  to  provide  possible  methods,  we  have  assumed,  not  quite 
correctly,  that  the  problem  of  designing  the  computation  procedure  may  be  split 
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into  two  or  three  separate  problems.  The  first  of  these  is  concerned  with  impor¬ 
tance  sampling,  taken  up  in  Section  2.  This  is  understood  to  mean,  for  our  pur¬ 
poses,  altering  the  natural  history  of  a  neutron  or  gamma  ray  so  as  to  insure  an 
adequate  supply  of  collisions  even  at  large  distances  from  the  source.  We  assume 
that  the  blast  wave  and  similar  complications  have  no  essential  role  in  designing 
the  importance  sampling. 

The  second  problem,  with  which  Section  3  is  concerned,  is  that  of  scoring. 
This  means  using  the  collisions  which  make  up  the  histories  to  get  answers  of  any 
desired  kind.  Considerable  emphasis  is  placed  upon  efficient  means  of  computing 
flux  or  dose  at  a  point.  Finally,  Section  4  is  devoted  to  a  discussion  of  methods 
of  eliminating  unncessary  flux  calculations  when  there  are  many  detectors. 

A'though  the  work  cannot  be  considered  complete  until  the  methods  are 
tried  out  in  the  final  code  or  a  similar  one,  specific  recommendations  are  made 
for  importance  sampling  and  scoring  in  Sections  2  and  3.  The  material  in  Sec¬ 
tion  4  has  not  been  completed  to  the  point  where  a  recommendation  can  be  made; 
however,  the  method  may  be  introduced  easily  into  any  code. 

Finally,  considerable  work  remains  in  reducing  the  methods  proposed  to 
practical  flow  charts  for  coding  into  the  transport  calculation. 
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2.  IMPORTANCE  SAMPLING 


2.1  GENERAL  THEORY 

In  order  to  discuss  the  aims  and  problems  of  importance  sampling,  let  us 
review  the  basic  theory.  Write  the  integral  Boltzmann  equation  for  x(P)dPi  the 
density  of  particles  coming  out  of  collision  in  dP  (in  phase  space)  as 

X  (P)  =  ;  K(P'^P)  X  (POdP  +  S(P).  (2.1) 

K(P'-P)  dP  is  the  e:q)ected  number  emerging  from  a  collision  at  P  in  dP,  given 
that  a  particle  came  out  of  collision  at  P'.  S(P)  is  the  source  density.  Eq.  2.1 
may  be  transformed  with  the  use  of  an  arbitrary  positive  function  I(P).  Multiply 
through  by  the  factor 

m 

!  I(P)  S(P)  dP 
and  define 

~  _  X(P)  I(P) 

^  /  I(P0  S(P')  dP'  ^  ^  ^ 


X(P)  =  /  K(P'-P) 


I(P) 

KPO 


x(P0  dP' 


^  S(P)J^ 

/  s(;^  i(p)  dP 


(2.3) 


This  equation  is  the  same,  formally,  as  Eq.  2.1,  with  a  new  normalized  source 
distribution,  and  a  new  kernel 


K(P'-P)  = 


m. 

I(P') 


K(P'-P) 


(2.4) 
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Suppose  that  a  single  answer  of  the  form 

<f>  =  /  x{P)  f(P)  dP  (2.5) 

is  required.  Then  it  can  be  shown®  that  a  particular  choice  of  J(P)  satisfying  the 
adjoint  equation 

J(P)  =  /  K(P-.P')  J(P0  dP'  +  f(P)  (2.6) 

used  for  I(P)  permits  the  answer  (Eq.  2.5)  to  be  obtained  exactly  (l.e.,  with  zero 
variance). 

The  function  J(P),  the  “importance  function,”  has  the  physical  meaning  of 
the  e;qpected  final  contribution  to  the  required  answer  <f>  for  a  particle  coming 
out  of  collision  at  P.  The  altered  kernel  has  the  reasonable  property  that  it 
should  be  biased  to  encourage  collisions  where  they  are  likely  to  give  significant 
contributions  to  the  answers.  Beyond  trying  to  exploit  this  last  feature  impor¬ 
tance  sampling  is  not  often  carried  out  because  of  the  following  difficulties. 

1 .  J(P)  is  never  Imown  exactly.  If  it  were  then  the  identity 

;  J(P)  S(P)  dP  =  /  f(P)  x(P)  dP  (2.7) 

would  permit  solution  by  quadrature.  Solution  of  the  adjoint  equation  is  as  diffi¬ 
cult  as  solution  of  the  original  equation.  However,  it  has  been  argued®  that  since 
the  variance  is  positive,  a  choice  I(P)  which  has  the  essential  features  of  J(P)  may 
make  the  variance  very  small.  The  physical  meaning  of  J(P)  can  often  be  used  to 
help  estimate  a  “reasonable”  guess  of  I(P). 

2.  K  we  attempt  to  use  an  I(P)  =  J(P)  the  normalization  of  the  kernel 
(Eq.  2.4)  is  unknown.  That  is 

N(P')  =  JdP  K(P'-P)  (2.8) 

can  not  be  evaluated  in  closed  form  and  the  expected  number  of  particles  emerg¬ 
ing  from  all  collisions  can  not  be  evaluated.  This  difficulty  can  be  evaded  in  two 
ways.  One  is  to  use  an  I  such  that  the  integral  (Eq.  2.8)  can  be  evaluated.  The 
exponential  transformation  may  be  regarded  as  such  a  choice.  Another  is  to 
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form  any  kernel  K'  (P'— P)  which  is  similar  to  K(P'— P),  pick  from  K',  evaluate 


N'  = 


k  (P'-P) 
K'(P'-P) 


and  create  this  number  of  particles  coming  out  of  collision  at  P.  This  generaliza¬ 
tion  of  the  usual  splitting  and  Russian  roulette  may,  in  turn,  be  carried  out  in  sev¬ 
eral  ways.  If  we  put  K'  =  K,  then  the  procedure  is  simply  to  track  particles  in  the 
usual  (unbiased)  way  from  collision  to  collision  and  create  I(P)/1(P0  particles  for 
each  one  that  comes  out  of  collision  at  P.  This  is  simple  and  direct,  but  it  has  the 
disadvantage  that  the  important  parts  of  phase  space  are  encountered  rarely  (if 
we  have  to  resort  to  importance  sampling  at  all).  Then  the  result  is  to  create 
many  particles  very  infrequently  -  a  procedure  that  is  characterized  by  a  large 
variance. 

A  somewhat  more  sophisticated  procedure  in  which  an  effort  was  made  to 
approximate  K  has  been  carried  out  successfully*  and  shown  to  permit  the  reliable 
computation  by  Monte  Carlo  methods  of  penetration  through  very  deep  layers. 

3.  In  many  problems,  there  is  no  single  answer  of  primary  interest.  On 
the  contrary,  a  map  of  the  radiation  field  as  a  function  of  position,  energy,  and 
time  may  be  desired.  If  importance  sampling  is  carried  out  as  outlined  above, 
one  particular  answer  (e.g.,  at  the  deepest  penetration)  may  be  estimated  well  at 
the  expense  of  accuracy  of  the  others. 

If  importance  sampling  could  be  designed  very  well,  then  it  may  happen  that 
a  set  of  different  estimates  of  the  same  radiation  field  may  be  better  obtained  by 
separate  computations  with  different  biasing.  It  is  likely,  for  instance,  that  sec¬ 
ondary  gamma  sources  in  the  atmosphere  are  best  derived  from  unbiased  neutron 
histories  separate  from  those  used  to  give  deep  neutron  penetration. 

On  the  other  hand,  if  we  assume  that  it  is  not  possible  to  design  the  impor¬ 
tance  sampling  so  that  individual  answers  are  given  very  precisely,  it  is  desirable 
to  introduce  an  importance  function  which  gives  reasonably  good  estimates  for 
detectors  spread  throughout  space.  We  suppose  that  these  detectors  are  located 
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In  the  vicinity  of  r^,  i  =  1,  2,  ,  and  that  the  total  flux,  integrated  over  time  and 

energy,  is  to  be  estimated  with  a  given  precision  at  each  one.  The  last  assumption 
does  not  preclude  the  possibility  of  recording  time  and  energy  dependence.  What 
it  means  is  that  where,  at  a  given  position,  the  flux  is  small  at  some  time  as  com¬ 
pared  with  a  different  time,  that  we  will  be  content  with  the  lower  precision  at¬ 
tached  to  the  estimate  of  the  smaller  flux. 

For  purposes  of  importance  sampling,  we  ask  that  the  several  fluxes  be  es¬ 
timated  in  such  a  way  that  the  weighted  sum 

<f>  =  S  Wi(f)(ri)  (2.9) 

is  estimated  with  minimum  variance.  If  we  wish  to  estimate  the  fluxes  so  that 
each  has  approximately  a  constant  relative  error,  then 

Wi  oc  [(f)  (£!)]■*.  (2.10) 

More  specifically,  we  should  make  wj  inversely  proportional  to  the  total  e^qjected 
contribution  from  which  a  flux  value  is  deduced.  For  example,  suppose  that  we 
wish  to  estimate  the  flux  in  a  number  of  spherical  shells  of  constant  thickness 
placed  about  a  point  source  of  radiation.  Then  for  each  shell  an  integral  propor¬ 
tional  to  4  7rr^  (f>  is  obtained,  the  prescription  for  weighting  is  then  that  we  should 
set 

1 

w  =  T - T-r 

4  7rr*<f) 

which,  in  turn,  may  be  estimated  roughly  as 
w  *  e^  . 

Since  we  wish  to  obtain  fluxes  weighted  by  w,  it  is  reasonable  to  use  w  as  an  im¬ 
portance  function,  I. 

It  is  worth  re-emphasizing  that  this  choice  of  importance  function  derives 
primarily  from  the  desire  to  estimate  fluxes  with  approximately  constant  relative 
error  throughout  a  certain  range,  rather  than  a  desire  to  obtain  penetration 
through  a  shield. 
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2.2  AN  EXPONENTIAL  TRANSFORMATION  FOR  SPHERICAL  GEOMETRY 

Let  us  Illustrate  the  general  notions  of  the  preceding  paragraph  with  a  some¬ 
what  simplified,  although  useful,  example.  Consider  the  transport  in  a  homogene¬ 
ous  infinite  medium  with  a  point  source.  The  spatial  dependence  of  the  flux  is 
through  the  radial  distance,  r,  alone. 

To  carry  out  importance  sampling  with  an  importance  function 

I  =  e^^  (2.11) 

we  split  the  kernel  K(P^-P)into  a  “transport”  and  a  “collision”  kernel; 

K(P'-P)  =  T(r'-r  I  E',nOC(E'-E,  n'-filr)  (2.12) 

where  T  gives  the  expected  number  of  particles  going  into  collision  in  dr  at  r 
given  that  a  particle  came  out  of  collision  at  r'  with  energy  E',  direction  Q'*  C  gives 
the  expected  number  of  particles  coming  out  of  collision  at  r  with  energy  in  dE, 
direction  in  dn,  given  that  one  went  into  collision  with  E',f2'.  With  an  importance 
function  that  depends  upon  position  only,  the  altered  kernel 

K(P'-P)  -  K(P'-P)  (2.13) 

requires  a  change  in  the  transport  kernel  only: 

k  =  TC  (2.14) 

T=I|^  T(r'«rlE',2)  .  (2.15) 

It  is  now  necessary  to  establish  a  random  walk  whose  transport  kernel  is  given 
by  T. 

It  has  been  pointed  out  that  this  is,  to  a  large  extent,  arbitrary.  One  method, 
for  example,  is  to  follow  particles  from  r'  to  r  using  the  kernel  T,  and  then  create, 
on  the  average,  I(r)/I(rO  particles  going  into  collision  at  r  .  This  method,  a  gen¬ 
eralization  of  splitting  and  Russian  roulette,  is  easy  to  carry  out.  It  has  the  im¬ 
portant  disadvantage  that  when  a  particle  flies  radially  outward,  then,  with  low 
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probability  many  particles  are  createtl.  Thus  additional  variance,  which  may  be 
large,  Is  associated  with  the  splltlng  process.  It  Is  desirable  to  use  a  kernel  as 
close  10  T  as  can  be  arranged  conveniently. 


For  the  spherical  problem  under  discussion,  the  appropriate  geometry  Is 
sketched  in  the  figure. 

Suppose  a  particle  starts  a  flight  at  B,  in  a  direction  whose  cosine  Is  o)  with 
respect  to  the  radius  from  the  source.  The  original  transport  kernel  is  given  by 
the  conditions  that 

a.  The  next  collision  is  along  the  ray  BC. 

b.  The  expected  number  of  collisions  in  ds  at  s  is  given  by 

T(s)  ds  =  n(E')  ds  (2.16) 

where  ii(E')  =  n*  is  the  macroscopic  cross  section  at  energy  E'. 

With  the  importance  function 

the  altered  transport  kernel  satisfies  (a)  above,  but  (b)  is  replaced  by 
f  (s)  =  T(s) 

=  H*  exp(-/s  +  fio  Vro  +  2r(,u)S  +  s^  -Mofo) 

=  /jL'exp[-A(s)]  (2.17) 

This  kernel  is  not  easy  to  pick  from  (for  one  thing  it  is  not  normalized  to  give 
exactly  one  next  collision).  It  is  straightforward,  however,  to  pick  from  the  kernel 
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T'(s)  =  +  exp  (-A(s)] . 

To  do  it,  pick  Xg  from  the  p.d.f.  e  ^  and  determine  Sg  from 
X(s#)=  Ag 

The  solution  for  Sg  is 

_  AgH^  -  Mo(4*-  tt>Mo)ro  +  a>^g)M■^  2(a)M*-  Mo)  Agrg  + Ag]^/^ 


Sg 


and 


r  =  (rg  +  2a)rgSg  +  Sg)'/^ 


Ag 

=  To  +  -  Sg  -  ^ 

Mo  Mo 


(2.18) 


(2.19) 


(2.20) 


In  doing  this  we  assume  that  ii*  >  Hq.  This  is  not,  of  course,  a  necessary  assump¬ 
tion.  To  relax  it,  however,  it  is  necessary  to  be  more  careful  about  the  behavior 
of  the  importance  function  for  large  r  and  to  take  care  to  calculate  Sg  correctly 
when  n'  =  Jig. 

Having  picked  s  from  a  kernel  other  than  T,  it  is  necessary  to  arrange  that 
N  particles  go  into  collision  where 


But 


N  =  — 

T'  ^ 

ds 


—  Mo(sH-raa>) 

Vr  J  +  2rgujs  +  s^ 


u^r 

N  - - - 

^  r  -  ^S(^-^u)Tq 


(2.21) 


-  Mo(s  +  rgw) 


=  +lx'  - 


(2.22) 


Since 


M'  >  Mo 


I 


and 


r  >  s  +  WFo 
we  have 

-4—  >  N  >  0  (2.23) 

The  factor  N  can  be  used  as  a  weight  factor  or  as  a  particle  multiplication  factor. 
The  latter  seems  more  natural.  In  a  straight  ahead,  one- velocity  transport  prob¬ 
lem  in  which  the  absorption  probability  is  1-C,  the  flux  behaves  like 

47rrV  =  exp  [-(1-C)  ^x'r] 

It  is  natural  to  take  =  (1-C)|i;  for  which 


Hence,  if  N  is  taken  to  be  a  particle  splitting  factor  as  discussed  here,  the  ex¬ 
pected  number  of  particles  coming  out  of  collision  is  NC  =1.  In  this  simple  prob¬ 
lem,  suppression  of  absorptions  comes  about  in  a  natural  way.  We  suppose  that  in 
a  real  spherical  calculation  that  the  same  behavior  will  hold  at  least  qualitatively 
and  that  N  should  be  used  to  create  extra  particles  rather  than  as  a  weight  factor. 
Naturally,  it  is  best  to  compute  N  collision  and  arrange  that  this 

many  particles  come  out  of  collision  there. 

The  e^qjonential  transformation  for  spherical  geometry  that  has  been  dis¬ 
cussed  here  may  be  generalized  to  concentric  spherical  shells.  It  will  not  be  dis¬ 
cussed  further,  since  it  has  been  included  primarily  for  purposes  of  illustration. 
The  geometrical  complication  associated  with  the  inhomogeneous  atmosphere 
makes  it  impossible  to  carry  out  importance  sampling  this  well.  Before  discussing 
the  atmosphere  question  in  detail,  we  carry  out  a  further  simplification  of  the 
spherical  problem  which  suggests  a  simple  treatment  of  the  atmosphere  impor¬ 
tance  sampling.  This  involves  making,  in  Eq.  2.17,  the  two  limiting  approxima¬ 
tions  for  small  and  large  s: 
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r  =  Vs^  +  iJojroS  +  rj 
»  To  +  ws  s  <  <  To 

s  s  +  wro  s>>ro  (2.24) 

These  are  exactly  right  when  w  =  ±1  and  otherwise  underestimate  r  by,  at  most, 
1/2  rj.  We  propose  to  use  these  approximations  in  T  to  obtain  a  T''. 

^  exp  [-Xi(s)]  (2.25) 


Xi(s)  =  ^.'s  -  iioTo  -  poojs  +  /ioro  s  rj 
=  n*s  -  MoS  -  s  ^  r 

To  pick  from  T'^,  let  \q  be  chosen  from  e 

and  if 


(2.26) 


Xo  (M'-a)po)ro; 


m'-wMo 


if 


Xo  >(^JL'-co^lo)  rj; 


_ Xo-fio(l~ci))ro 

S  “  . 

M-Po 


(2.27) 


Effectively,  these  correspond  to  carrying  out  the  exponential  transformation 
for  the  plane  for  s  <  rg,  and  assuming  a  straight  ahead  approximation  for  s  >  rQ. 

N  is  now  computed  as 


p'-ojpo 


exp  [(Jo(r-ro)-Pows]; 


s  ^  ro. 


P^-Pfl 


exp  [po(r-s-u)ro)]; 


s  >  ro 


(2.28) 


r  =  (ro  + 2a)roS +8^)^^^ 

Wells^  has  suggested  using  the  first  part  of  Eq.  2.27  alone  in  spherical 
problems. 
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2.3  IMPORTANCE  SAMPUNG  IN  AN  INHOMOGENEOUS  ATMOSPHERE 

The  change  in  atmospheric  density  with  altitude  complicates  any  attempt  at 
Importance  sampling  in  two  ways;  as  we  have  seen,  it  is  necessary  to  introduce  a 
guess  at  the  nature  of  the  solution  and  this  is  more  difficult.  Carrying  out  the  im¬ 
plied  importance  sampling  is  also  more  involved. 

A  relatively  straightforward  method  of  importance  sampling  follows  if  we 
subdivide  space  into  a  large  number  of  finite  volumes.  For  an  atmosphere  in 
which  the  density  depends  only  upon  vertical  height,  these  volumes  may  be  defined 
by  the  intersection  of  horizontal  planes  and  concentric  cylinders  with  a  vertical 
axis.  We  assume  an  importance  function  which  is  piece- wise  constant  in  each  of 
these  volumes.  The  values  may  be  computed  to  be  inversely  proportional  to  an 
integrated  flux  as  indicated  in  the  preceding  section.  On  the  other  hand,  the  values 
may  be  chosen  so  as  to  emphasize  regions  of  space  close  to  particular  detectors. 

In  any  case,  this  importance  function  may  be  used  to  define  a  splitting-and- 
Russian- roulette  procedure:  when  a  particle  is  tracked  from  a  volume  in  which 
the  importance  function  is  Ii  to  a  volume  in  which  it  is  I2,  arrange  that 

N  -  I2/I1 

particles,  on  the  average  actually  enter  the  new  region.  If  N  <  1,  N  is  interpreted 
as  a  probability  of  continuing  the  history  (1-N  is  the  probability  of  Russian  rou¬ 
lette).  If  N  >  1,  compute  the  integer  part  {N}  and  fractional  part  N-{N}  respective¬ 
ly.  Arrange  that  {n}  plus  with  probability  N-{n}  one  more  particles  actually  enter 
the  region.  This  is  usually  done  by  following  one  particle,  and  enterii^  coordi¬ 
nates  of  the  remaining  particles  in  a  ‘‘latent”  list.  The  particles  in  a  “latent” 
list  are  taken  up  one  at  a  time  after  the  preceding  parts  of  a  history  are  termi¬ 
nated  in  the  usual  way. 

This  splitting  provides  a  very  flexible  and  straightforward  way  of  increas¬ 
ing  particle  populations  in  otherwise  unlikely  parts  of  space.  It  has  been  used 
successfully  in  a  number  of  applications. 
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The  disadvantages  of  the  method  are  as  follows.  First,  for  the  atmosphere 
problem,  it  requires  the  introduction  of  an  otherwise  extraneous  geometry  with  the 
necessity  of  tracking  through  it.  Second,  it  requires  the  specification  of  the  many 
values  of  the  importance  function.  Third,  the  process  of  splitting  introduces  addi¬ 
tional  variance.  The  third  objection  should  not  be  understood  to  mean  that  splitting 
is  not  worthwhile,  but  that  alternative  ways  of  using  the  importance  function  may 
be  better.*  The  second  objection  may  be  overcome  if  an  analytic  computation  of 
the  importance  function  can  be  devised. 

In  view  of  these  disadvantages,  especially  the  first,  we  have  given  some  at¬ 
tention  to  the  possibility  of  devising  an  alternative  scheme,  more  nearly  along  the 
lines  of  the  earlier  section.  The  first  part  of  such  a  scheme  requires  the  specifi¬ 
cation  of  a  simple  approximate  form  for  the  importance  function. 

A  simple  generalization  of  the  radial  exponential  is 


I(r)  =  exp  [Ofl  /  p(s)  ds] 

(2.29) 

where  the  integral  of  the  density  is  taken  along  the  line  from  source  to  point  r. 

introduce  the  cumulative  density  function 

F(z)  =  r  p(z)  dz 

J  z 

(2.30) 

where  z  is  vertical  height  and  p(z)  is  the  atmospheric  density  at  z. 

specified  by  (x,y,z)  and 

Let  r  be 

y  =  z/r  =  z(x^  +  y^+zV''^^ 

(2.31) 

I(x,y,z)  =  exp|^  [F(zo)-F(z)]| 

(2.32) 

*One  such  method  has  been  proposed:  it  consists  in  examining  the  segments  cut 
off  the  particle  path  by  the  volume  subdivision,  and  computing  the  expected  num¬ 
ber  of  particles  which  should  go  into  collision  in  each  segment.  This  number  is 
given  by  (I  /  pe'^ds)  for  each  segment.  One  arranges  that  this  number,  on  the 
average,  go  into  collision  in  that  segment.  The  method  has  the  apparent  advan¬ 
tage  that  long  flights  do  not  depend  upon  the  survival  of  a  cascade  through  the 
shorter  distances.  It  has  not  been  tried  to  any  extent. 
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where  the  source  is  at  (0,0, z^).  The  value  of  Oq  is  adjustable. 

This  form  of  the  importance  function  must  be  regarded  as  tentative.  It  is 
reasonable  when  the  variation  in  density  is  not  large  since  it  reduces  to  exp  (-a^r) 
in  that  case.  For  very  large  distances  when  the  density  is  rapidly  varying  it  may 
well  be  wrong.  For  example,  it  may  be  easiest  to  achieve  deep  penetration  by  a 
flight  to  a  thin  part  of  the  atmosphere,  followed  by  a  long  flight  through  this  region. 
An  “indirect”  penetration  of  this  kind  is  not  properly  emphasized  by  the  impor¬ 
tance  function. 

The  validity  of  the  proposed  function  has  been  tested  in  two  ways.  First, 
very  crude  multiple  scattering  calculations  were  made  to  see  whether  an  indirect 
flight  can  lead  to  much  larger  fluxes  at  great  distances  than  a  direct  flight.  In 
these  calculations  distances  were  chosen  so  that  0.1  rad  would  result  from  a 
100  megaton  yield.  This  pushes  the  parameters  to  their  limit.  An  effective  atten¬ 
uation  coefficient  of  0.04  cmVg  was  used  throughout.  Multiple  scattering  in  vol¬ 
umes  about  one  mean  free  path  on  a  side  was  estimated  for  various  positions  of 
source  and  scattering  volumes.  According  to  the  rough  estimates  obtained,  the 
indirect  scattering  can  exceed  the  direct  flux  when  the  source  is  above  12  km. 

Even  at  the  extreme  ranges,  the  excess  is,  at  most,  a  factor  of  100.  This  ratio  is, 
of  course,  very  important  itself,  but  is  not  considered  to  be  a  serious  error  in  de¬ 
signing  an  importance  function. 

A  second  check  on  the  flux  estimate  was  made  by  plotting  4nr^(f>  vs  /pds 
from  some  calculations  already  made  of  penetrations  in  the  atmosphere.  These 
calculations  may  not  be  entirely  reliable  at  very  deep  penetrations,  since  no  im¬ 
portance  sampling  was  involved. 

In  Figs.  2.1  and  2.2  the  results  for  total  flux  in  the  highest  energy  bins  are 
plotted  for  two  source  energies,  two  source  heights,  and  several  detector  posi¬ 
tions.  The  points  do  not  fall  far  from  straight  lines,  indicating  again  that  the  flux 
form  in  Eq.  2.29  is  reasonable.  The  reason  for  the  difference  in  slopes  between 
the  “DOFL”  and  “SANDIA”  curves  in  Fig.  2.2  in  unknown,  but  is  presumed  to 
result  from  different  cross  sections. 


14 


Integrated  Density,  /pds,froni  Source  to  Detector  in  g/cm* 


Fig.  2.1  —  Neutron  penetration  in  the  atmosphere;  flux  of  1.5  to 
2  Mev  neutrons  from  a  2-Mev  source. 


IS 


Integrated  Density,  /pds,  from  Source  to  Detector,  in  g/cm* 


Fig.  2,2  —  Neutron  penetration  in  the  atmosphere;  flux  of  10  to 
14  Mev  neutrons  from  a  14-Mev  source. 
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Taking  Eq.  2.32  for  the  importance  function,  an  altered  transport  kernel  Is 


defined  as  before 

T(e'-E)=-|s|  T(r'-r) 

(2.33) 

^  (F(z.)  -F(z)J  -  ^  (F(zJ  -  F(z')]| 

(2.34) 

where 

t 

’  F- 

Aside  from  necessary  delta  functions, 

T(r'-r  lEOds  =  a(E0p(z)  expj-  ^^^[F(zO-F(z)]|ds  (2.35) 

where  s  =  l£“£^  I  is  the  distance  between  r  and  r'  as  the  neutron  or  gamma  flies. 
Then 

T  =  a(E0p(z)  exp|-^|^^[F(z')-F(z)]+^  [F(zo)-  F(z)]-^  [F(zo)-  F(zO]l 

(2.36) 

where  ujz  =  (z-z^/s. 

We  have  been  unable  to  devise  a  T'  which  resembles  T  as  closely  as  was 
done  before  and  which  permits  an  easy  determination  of  the  distance  to  be  trav¬ 
eled.  If,  as  before,  we  set  the  exponent  in  T  equal  to  -Ao  (chosen  from  e  ~^)  we  get 
a  transcendental  equation  for  z: 

Afl  =  -  ^  F(z)  +  —  F(z)  +  const  (2.37) 

Y 

Note  that  a(E')/co2  ^  constant  which  is  known  but  y  depends  upon  z,  and  that 

F(z)  is  approximately  exp(-Az).  If  the  last  equation  is  easily  solved  (e.g.,  by  trial 
and  error),  and  if  dA/ds  *  0,  then  a  scheme  like  the  exponential  transformation 
can  be  worked  out.  We  have  assumed  that  this  would  be  inconvenient  and  propose 
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a  direct  ad  hoc  application  of  the  approximate  exponential  transformations  of 
Eq.  2.27.  That  is,  define  a  T''  using 


a(E')  for  s  “S  r' 

and 

a(E')  -  Co  for  s  >  r'; 

here 

r'‘(r-r') 

O)  = - ; - 

r's 

so  that 

T''(s)  ds  =  p  (z)[a(E0-  wuojds  exp  [  W)  -  F(z)]|; 

=  p(z)  [a(EO-Uo]  ds  exp  •[-  [F(zi)-F(z)] 

I 

+  [F(zo)-F(zi)]  I ;  s  >  r' 

Zi  is  that  value  of  z  for  which  s  =  r\  i.e., 

Zj  =  z'  +  r' 

To  carry  this  out,  pick  a  X  from  e  ^  and 

1.  If  X  ^  =  [F(zo-F(zi)]  , 

<^z 

z  is  determined  from 
F(z)  =  F(z,)  - 

and 


N  = 


a'-<j}(XQ 


exp{-c.(i-  i)[F(z.)  -  F(z')]-  5,^ 


+ ; 


(2.38) 


s  <  r'. 


(2.39) 


(2.40) 


(2.41) 

}■ 

(2.42) 


II 


2.  If  A,  >  Xj,  z  is  determined  from 


F(z)  =  F(z,)  - 


u>2  (X-A^) 
a'-CTo 


and 


N  = 


exp 


(i  -  i)  [FW-  F(z.,l-(^  - 


(2.43) 

It  is  recommended  that  if  N  particles  are  to  go  into  the  next  collision,  and 
if  ag/aT  is  the  probability  of  survival,  that  it  be  arranged  that  just  N  os/crx>  the 
average,  come  out  of  the  next  collision.  Naturally,  the  proper  proportion  of  elastic, 
inelastic  scatterings,  and  other  possible  interactions  must  be  preserved. 

Finally,  some  means  must  be  employed  to  cut  off  the  importance  sampling  at 
very  large  distances  so  as  to  prevent  particles  from  wandering  out  to  infinity.  A 
simple  method  of  doing  this  consists  in  terminating  a  history  when  a  particle 
leaves  a  large  volume  which  encloses  the  problem.  The  cutoff  should  be  suffi¬ 
ciently  far  from  a  detector  that  it  has  a  negligible  effect  on  the  flux  at  that  detector. 
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3.  SCORING  METHODS  -  FLUX  AT  A  POINT 


3.1  UNCOLLIDED  FLUX  ESTIMATOR 

In  some  calculations  of  a  radiation  field,  it  is  desirable  to  obtain  answers 
for  flux,  dose,  or  the  like  at  a  point.  The  alternative  -  to  estimate  similar  quan¬ 
tities  as  averaged  over  finite  volumes  -  may  be  satisfactory  if  some  care  is  used 
in  choosing  the  volume  small  enough  so  that  the  difference  between  averaged  an¬ 
swer  and  its  point  value  is  small.  If  this  volume  is  very  small,  then  there  is  no 
essential  difference  from  flux  at  a  point. 

Since  no  history  can  be  expected  to  carry  a  neutron  through  a  point,  it  is 
necessary  to  use  analytic  expected  value  estimates.  One  common  method  uses 
the  relation  between  il(r\E,Q)  dr'dE  d^,  the  density  of  particles  entering  collision 
in  volume  dr  with  energy  in  dE  and  direction  in  dfl,  and  the  flux  at  r: 


f  r—r^ 

(pir)  =  /  dE'  /  dr'  f  dE  f  dQ  4,  (r',E,Q)  x  ^  2g  f  E 


^-fi(E')  Ir-r' 
4n  Ir-r'l*^ 


In  this  expression  g(aj,E)  dw  is  the  probability  of  scattering  a  neutron  of  energy  E 
through  an  angle  whose  cosine  is  in  the  range  (w,  w+doj).  E'  is  the  energy  after 
scattering  and  m(E')  is  the  total  attenuation  coefficient  of  the  scattered  neutron. 
Isotropic  scattering  gives 


_  1  OS 

®  2  O'j' 


(3.2) 
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The  straightforward  method  of  using  this  expression  is  to  evaluate  the  esti¬ 
mator  in  braces  every  time  a  particle  goes  into  collision  and  average  over  his¬ 
tories. 

It  is  easy  to  see  that  this  procedure  yields  an  infinite  variance  if  the  point  r 
lies  within  the  medium.  In  the  neighborhood  of  r',  let  r  =  Ir-r^  I;  the  integral  is 
essentially  of  the  form 

^  ~  ^(rO  /  rW  (3.3) 

which  converges.  The  mean  square  flux  has  contributions  like 

^  ~  4){r*)  j  rMr  (3.4) 

which  does  not  converge. 

The  fact  that  the  variance  is  infinite  does  not  necessarily  imply  that  the  use 
of  this  procedure  is  improper.  Indeed  it  will  be  proved  that  for  many  histories, 
the  mean  value  of  the  flux  obtained  in  this  way  converges  to  the  correct  answer. 

On  the  other  hand,  the  analysis  shows  that  the  infinite  variance  method  has  an  er¬ 
ror  that  converges  more  slowly  than  if  the  variance  is  finite.  For  this  reason  it 
is  desirable  to  find  a  method  which  has  a  finite  and  reasonably  small  variance. 

This  section  will  be  devoted  to  an  analysis  of  the  effects  of  the  infinite  variance 
and  a  discussion  of  several  methods  of  making  the  variance  finite. 

3.2  ANALYSIS  OF  THE  INFINITE  VARIANCE  ESTIMATOR 

According  to  the  previous  analysis,  the  infinity  in  the  variance  of  the  un¬ 
collided  flux  estimator  arises  from  the  behavior  of  the  estimator  near  r  =  0, 
namely  as  r"^  Essentially  the  same  situation  should  occur  if  we  sample  a  random 
variable  x  with  a  p.d.f.  oc  near  the  origin  and  evaluate  an  estimator  proportional 
to  x"^.  Let  us  analyze  this  situation  where,  for  convenience,  the  range  of  x  is 
(0,1). 
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Suppose 

p(x)  =  3x^  0  ^  X  <  1 


(3.5) 


and  the  estimator  is 
z(x)  =  (3xV* 


Clearly 

<z> 


<  z^> 


p(x)  z(x)  dx  =  1 


p(x)  [z(x)]2  dx  = 


OO 


(3.6) 


(3.7) 


To  determine  the  convergence  to  the  mean  of  the  average  value  of  z,  let  us  first 
find  the  p.d.f.  of  z: 


Pz(z)  =  - 


dx 


dz 


27T 


1 

-^Z<oo. 


(3.8) 


In  order  to  calculate  the  behavior  of  the  sum  of  many  values  of  z  chosen 
independently,  a  knowledge  of  the  characteristic  function*  is  required.  The  latter 
is 


ch2(t)  — 


e‘^‘^Pz(z)  dz  = 


z”*/^  e^^^dz 


(3.9) 


This  integral  is  not  easy  to  determine  in  closed  form.  However,  its  behavior 
near  the  origin  may  be  computed.  Let 

h(t)  =  r*/2[chz(t) -  1  -  it]  (3.10) 


itz 

e 


-I  -itz 


dz 
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(3.11) 


1 

2/3" 


eiu  -i-iu 
■ 


e^“-l-iu, 
- du 


Now,  since  the  integrand  behaves  like  u"^/^  near  the  origin,  the  first  integral  ex¬ 
ists  and  the  second  is  proportional  to  t*/^.  Thus,  near  t  =  0, 

chz(t)  =  1  +  it  +  ytV^  +  St^  +  .  .  .  (3.12) 

log  chz(t)  =  it  +  ytV2  +  S't^  +  .  .  . 


Suppose  a  z  is  computed  independently  n  times,  and 


n 


Z  =  2  Zi 
i=l 


(3.13) 


n 

chz(t)  =  n  =  [chz(t)]"  (3.14) 

i=l 

log  chz(t)  =  n  log  ch2(t)  =  int  +  ynt®/^  +  6'nt^  +  .  .  .  (3.15) 

When  a  linear  change  is  made  in  a  random  variable: 

?  =  a  +  bZ  (3.16) 

we  get 

ch^(t)  ^  /  e‘*^-  pj(?)  d? 

=  /  eit(a+bZ)p^(Z)dZ 

=  eitachz(bt)  (3.17) 

Then  if  we  put 

?  =  ^^  =  -n‘/'  +  n-2/3z  (3.18) 

log  ch^(t)  =  -in‘/®t  +  in  (n"VH)  +  yn  (n"^/®t)^/^  +  6*n(n"V®t)^  +  .  .  . 

-  yt®/2  +  6'n-‘/® t^  +  .  .  .  (3.19) 
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In  the  limit  n  —  «> , 


log  chf(t)  =  yt®/^  (3.20) 

The  inversion  of  this  characteristic  function  leads  to  a  unique  p.d.f.  for  £.  At  the 
moment,  the  precise  form  of  this  is  not  important.  Note,  however,  that  the  relation 
defining  ?  in  terms  of  Z  has  the  form  of 

(3.21) 

where "z  denotes  an  “experimental”  —  i.e.,  Monte  Carlo  -  mean.  This  means  that 
the  difference"?  -  <z>  has  expected  value  zero,  its  deviation  from  zero  is  random¬ 
ly  distributed  according  to  a  p.d.f.  whose  width,  as  n  -  <»,  is  measured  by  n"‘/®. 
Finaliy  this  p.d.f.  behaves,  for  large  f,  as 

From  these  results,  it  may  be  concluded  that  in  this  situation,  the  fact  that 
the  variance  is  infinite  is  in  itself  no  reason  against  its  use.  The  first  require¬ 
ment  on  any  method,  that  after  repeated  sampling,  it  give  an  answer  close  to  the 
correct  answer,  is  preserved.  On  the  other  hand,  its  convergence  is  slower  than 
if  the  variance  exists  and  large  departures  from  the  expected  deviation  from  the 
mean  may  occur. 

It  may  be  objected  that  in  any  actual  calculation,  it  will  be  necessary  to 
truncate  z  and  that  this,  in  fact,  guarantees  a  finite  variance.  To  investigate  this, 
let  us  suppose  that  the  c.d.f.  of  z  is  truncated  in  the  following  wa3r. 


Ptz(z) 


1 

2/3’ 


Z/-5/2  =  1  - 


1 

3/3" 


1 

3  ^  z  <  Zo, 


Ptz(z)  =  1  - 


1 

3/3" 


Zo  «  z  <  3zo, 


Ptz(^)  =  1;  3zo  <  z  <  Zq. 

The  step  at  3z0  is  arranged  so  that 


(3.22) 
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I  zdPtz(z)  =  1 

*^i/Z 

Then  it  is  easy  to  see  that  <z^>  oc  <z*>  oc  Zq®'^^.  By  computing  the 

characteristic  function  of  z,  and  of  Z  =  E  zj  and  finally  of 

ft  =  ^-nV2 

it  is  possible  to  show  that,  although  the  variance  is  finite  (a  ZflV^),  it  requires  a 
sample  size  at  least  as  large  as  n  >  before  the  distribution  of  the  mean  be¬ 
comes  normal.  This  supports  the  intuitive  expectation  that  truncation  does  not 
help  as  long  as  the  expected  number  of  samples  in  the  truncated  part  of  the  distri¬ 
bution  is  small.  In  other  words,  if  the  truncation  is  to  do  any  good,  it  must  have  an 
appreciable  effect  on  the  average  result. 

To  return  to  the  estimation  of  flux  at  a  point,  we  note  that  the  truncation  can 
be  carried  out.  At  best,  however,  it  results  in  estimation  of  average  flux  over  a 
small  region.  The  argument  above  shows  that  this  region  should  be  chosen  so 
that  an  appreciable  number  of  collisions  take  place  inside  it  during  the  actual 
sampling  if  the  benefits  of  finite  variance  are  to  be  felt.  This,  of  course,  leads  to 
the  same  dilemma  in  choosing  volumes  as  one  has  in  the  first  place. 

3.3  FINITE  VARIANCE  METHODS 

As  a  general  rule,  to  reduce  the  variance  of  the  Monte  Carlo  estimate  of  an 
integral  it  is  desirable  to  sample  from  a  density  function  which  resembles,  as 
closely  as  possible,  the  entire  integrand.  If  the  integrand  includes  any  singular 
part,  it  should  be  incorporated  in  the  sampling  p.d.f.,  not  in  the  scoring  function. 

In  the  flux  estimation,  the  singularity  (in  the  second  moment)  derives  from 
the  r"^  behavior  of  the  estimator.  If  this  is  to  be  removed,  the  p.d.f,  which  deter¬ 
mines  the  point  at  which  the  collision  occurs  (the  collision  from  which  the  flux 
estimate  is  made),  must  be  altered  so  as  to  include  explicitly  the  r"^.  This  means 
in  turn  that  these  last  collisions  before  estimation  cannot  occur  in  the  usual  way. 
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What  is  required  is  that  for  every  collision,  we  estimate  the  “once-more  collided” 
contribution  to  the  flux  at  a  point,  and  choose  the  position  of  the  intermediate  col¬ 
lision  so  as  to  make  the  variance  finite.  If  this  procedure  is  carried  out  for  every 
collision  -  including  the  source  —  then  the  final  estimate  is  for  collided  flux  only. 

A  separate  analytic  calculation  must  be  made  of  the  uncollided  flux. 

To  determine  a  suitable  p.d.f.  for  the  intermediate  collision  point,  we  con¬ 
sider  the  integral  for  the  once  scattered  flux  at  a  detector  a  distance  R  away  from 
a  collision.  Assume,  for  simplicity,  that  scattering  is  isotropic,  that  the  cross 
section  is  unchanged  by  collision  and  that  the  scattering  probability  is  C. 


The  sketch  shows  the  geometry  with  the  collision  at  point  C,  the  detector 
at  D.  The  once  scattered  flux  at  D,  given  that  a  particle  goes  into*  collision  at  C 
is 


(f)i=  ;  dV 


-pTi 

e 

4Trr'f 


(3.23) 


The  singularity  in  the  variance  may  be  removed  if  the  p.d.f.  of  the  inter¬ 
mediate  point  A  contains  the  factor  dV/r^  rl .  Both  radii  are  needed;  if  only  r^  is 
used,  then  the  collision  points  occurring  close  to  C  would  again  give  infinite  vari¬ 
ance.  The  simplest  such  distribution  is 


Ui(R)  dV-p 


R  dV 


rrr 


1  ^2 


(3.24) 


*Note  that  if  a  particle  is  known  to  come  out  of  a  collision  at  C,  as  for  example 
in  the  source,  the  factor  in  this  integral  is  replaced  by  C. 
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If  the  intermediate  point  is  sampled  from  Ui(r),  then  the  flux  estimator  is 


C2  ^-n(r.*r,) 

Ui  (47r)'‘rtr|  16  R 


(3.25) 


We  note  that  when  C  falls  close  to  D,  R  becomes  small  and  fj  large.  How¬ 
ever,  since  the  singularity  is  only  R"*,  the  mean  square  of  fj  exists.  Higher  mo¬ 
ments  of  fi  diverge.  This  fact  does  not  negate  the  usual  form  of  the  central  limit 
theorem  but  it  means  that  Monte  Carlo  estimates  of  the  variance  itself  are  more 
than  usually  uncertain.  Moreover,  the  convergence  to  the  normal  law  will  be 
slower  than  if  <fi>  existed.  Higher  moments  of  the  scattering  may  be  made  to 
converge  by  including  more  orders  of  scattering  in  the  estimation  process. 

Although  Ui  removes  the  difficulty  associated  with  small  rj  or  r2,  its  be¬ 
havior  for  large  rj  and  r2  is  not  very  close  to  the  behavior  of  the  integrand  of  the 
single  scattering  integral.  The  following*  p.d.f.  may  be  better  in  that  respect: 


U2dV  =  ^ 
^  8w 


•  ■"■5 —  + 


‘1 


dV 


(3.26) 


For  fixed  R  and  r2  small. 


U2 


8tt 


2  2 

Ii+de  ^ 


ri  r2 


ifl  ^ 

8rr  r| 


(3.27) 


This  shows  that,  again  for  fixed  R,  the  use  of  U2  makes  the  single  scattering  inte¬ 
gration  have  finite  variance. 

For  R  small  the  mean  square  may  be  shown  to  vary  as 


ST!  =  (3.28) 

As  we  integrate  over  R  with  weight  factor  R^  the  final  integral  diverges  unless 
p'  is  made  to  vary  as  1/R  for  small  R.  Up  to  now  li'  is  arbitrary  so  that  it  may 
be  made  a  function  of  R  such  that 


*We  are  limiting  the  choice  of  U  to  functions  for  which  the  choice  is  straightfor¬ 
ward. 
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lim  R/i'(R) 
R-0 


(3.29) 


exists.  For  example,  we  may  use 


/  — 


(3.30) 


3.4  SAMPLING  FROM  Ui(r)  AND  Uzir) 
To  pick  a  point  from 


U.(r)  dV  =  ^ 


dV 


A 


use  u>  and  r]  as  coordinates; 
u)  is  the  cosine  of  the  angle 
between  rj  and  R. 


Then  write 


Ui(ri,a;)  dV- 


Ui(ri,c«;)  da;  d(p 


27r 


dc/^Uj  (rj.a;) 


XI  ^  2ti 

do;  I  dc/;Ui(ri,a))rfdri 
1  •'0 


The  integral 

M(ri)  =  r?  /  da;  /  dp  Ui(r„a;) 


(3.31) 

(3.32) 


is  the  marginal  distribution  for  rj.  If  we  pick  an  rj  from  M(ri)  drj,  then  the 
quantity  in  brackets  is  a  conditional  distribution  for  a;  and  (p  given  r^. 

A  straightforward  integration  shows  that 


2  1 

M(ri)  =zz— 

Ti  r, 


ri+R 

it;^ 


(3.33) 


28 


If  we  set 


ri  =  Riy. 

then  the  p.d.f.  of  y  is 

My(y)=pi  0^y< 


To  pick  from  this  function,  note  that 

*1/Y 

ry  OQ 


Jr.«.  ^1/Y 

My(y)  dy  =  I  My(yO  dy' 
Y 


and,  in  particuiar, 


Jmaa  y»l 

My(y)  dy  =  I  My(y)  dy 
1  •'0 


(3.34) 


(3.35) 


(3.36) 


Together  the  last  two  equations  show  that  to  pick  y  on  the  range  (0,00)  we  may 
pick  a  y'  from  My(yO  on  (0,1)  and  with  probabiiity  1/2  set  y  =  y',  otherwise 

y  =  l/y^ 

To  pick,  in  turn,  on  the  range  (0,1),  the  function  Y(^)  defined  by 
2  J  My(y)  dy  =  ^  (3.37) 


was  tabulated  and  fitted  by  a  rational  fraction: 
1.23  3  7  4  --  0.412593 


Y(0- 


1  -  0,17889 


The  maximum  relative  error  is  0.002,  which  seems  adequate. 
When  y  is  found,  co  must  be  chosen  from  the  p.d.f. 

,2 


w(„)  =  =  — y 


1 


M, 


log 


1+y  l+y^-2yci> 


ll-yl 


(3.38) 


(3.39) 
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The  solution  of 


(3.40) 


(3.41) 

(3.42) 

If  necessary,  an  azimuthal  angle  may  be  chosen  (at  random)  so  as  to  specify 
A,  the  point  of  intermediate  collision,  completely. 

To  pick  from  U2,  note  that  U2  is  an  average  of  two  p.d.f.’s  each  of  which  has 
the  form 


4n-  r^ 


Picking  from  U2  consists  simply  of  the  following:  (a)  pick  (with  equal  proba¬ 
bility)  either  the  collision  point  or  the  detector  point  as  an  origin;  (b)  pick  a  direc¬ 
tion  isotropically;  and  (c)  pick  a  distance  r  along  that  direction  from  the  p.d.f. 

-jLL^r 

li'e  .  This  distance  is  either  or  r2  and  the  remaining  distance  must  be  com¬ 
puted. 

The  advantages  in  using  U2  are  that  it  is  easy  to  sample,  and  that  it  is  possi¬ 
ble  to  generalize  in  a  straightforward  way.  We  may  use  various  linear  combina¬ 
tions  of  a  function  of  (ri,a)i)  and  a  function  of  (r2,a)2).  Including  u)  dependence 
means  that  the  direction  from  either  origin  is  chosen  in  a  biased  fashion. 


log 


1+y 

ll-yl 


do)'  =  I 


IS 


_  i+y^-(i+y)^[(i-y)V(i+y)^]^ 

2y 


and 


r2  =  VrffllJs-Zrjlla)  =  R(l+y) 


1-y 


1+y 
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A 


Note  that  if  the  once  collided  estimator,  fi(ri,r2)  is  a  symmetric  function  of 
rj  and  r^  (as  it  is  for  the  one  velocity  isotropic  scattering  case  described  above), 
the  following  simplification  occurs: 


/■ 


dV  U2(ri,r2)  f(ri,r2) 


Av 

87T  J 


fi(ri,r2) 


(3.43) 
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^2 


-Il'Vi 

^1 


fi(r2.ri) 

(3.44) 


This  means  that  the  choice  of  points  C  or  D  as  an  origin  may  be  dispensed  with. 

3.5  TESTS  OF  FLUX  ESTIMATORS 

A  series  of  computer  experiments  was  run  in  an  attempt  to  compare  the  re¬ 
sults  of  the  different  methods  of  scoring.  The  problem  solved  was  the  collided 
flux  distribution,  6^,  in  an  infinite  homogeneous  medium  which  scatters  isotropic¬ 
ally.  For  a  point  isotropic  source,  accurate  numerical  solutions  are  known. ^  A 
code  was  set  up  to  generate  histories  in  a  completely  straightforward  way.  That 
is,  no  importance  sampling  or  other  biasing  was  introduced.  Four  alternative 
scoring  methods  were  used: 

e-^r 

1.  Uncollided  flux  estimator:  -r—r 

47rr 

2a.  Once  collided  flux  estimator,  using  Uj 
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2b.  Once  collided  flux  estimator,  using  Uj 

3.  Estimating  flux  from  the  average  track  length  per  unit  volume  in 
a  thin  shell. 

The  answers  for  the  collided  flux  are  shown  in  Table  3.1,  along  with  the 
following  additional  data:  the  time  per  history,  the  number  of  histories,  the  ex¬ 
perimental  estimate  of  the  standard  deviation.  The  more  elaborate  methods,  2a 
and  2b,  take  more  time  per  history,  but  seem  to  have  smaller  standard  deviations 
for  a  fixed  number  of  histories.  A  criterion  which  combines  the  two  effects  is  de¬ 
rived  from  the  quantity 

Q  =  T  X  (variance),  (3,45) 

the  product  of  the  total  time  by  the  square  of  the  standard  deviation.  This  param¬ 
eter  should  be,  on  the  average,  independent  of  the  number  of  histories,  and  show 
whether  a  more  time-consuming  method  pays  off  in  more  accurate  answers.  This 
Q  is  also  tabulated  for  methods  2  and  3.  In  method  1,  the  variance  is  infinite  and 
Q  has  no  meaning. 

It  requires  considerable  sampling  to  get  a  variance  estimate  that  is  accurate 
enough  to  draw  good  curves  of  Q  as  a  function  of  distance  and  for  different  sam¬ 
pling  methods.  The  results  obtained  were  not  generally  as  accurate  as  might  be 
desired.  Nevertheless  we  believe  that  the  trends  shown  are  correct.  Method  2a 
seems  best,  with  2b  somewhat,  although  not  much,  worse.  Method  3  (track  length 
in  a  thin  shell)  generally  has  the  largest  Q.  For  C  =  0.9  the  differences  are  not 
great  and  in  two  problems,  method  3  was  better.  It  should  be  recalled,  however, 
that  in  method  3,  an  entire  shell  was  used  to  obtain  answers.  That  is,  full  use  was 
made  of  the  symmetry  of  the  problem.  This  is  entirely  justified  in  the  test  prob¬ 
lem,  but  not  in  any  calculation  without  the  symmetry.  Substantially  more  time 
would  be  required  to  get  answers  to  this  accuracy  in  a  small  area  on  the  surface 
of  the  shell  so  as  to  approximate  a  point.  To  be  more  specific,  suppose  that  a 
sampling  area  of  1  (mean  free  path)^  on  the  surface  is  used.  Then  it  would  take 
47rr*  times  as  long,  without  biasing,  to  get  equal  statistical  precision  as  compared 
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Table  3,1  —  Fiux  at  a  Point  -  A  Compaxison  of  Several  Methods 
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with  the  situation  where  the  whole  shell  is  used.  This  increases  Q  by  a  factor  of 
4ffr^  The  result  is  shown  in  the  last  column.  This  Q  is  much  iarger  than  the 
others. 

In  judging  the  usefulness  of  alternative  scoring  schemes  in  other  codes, 
some  care  must  be  used  in  interpreting  these  results.  In  the  tests,  the  time  re¬ 
quired  to  generate  a  singie  collision  was  reiatively  short.  In  a  code  in  which  a 
larger  fraction  of  time  is  devoted  to  tracking,  the  more  elaborate  and  slower 
methods,  2a  and  2b,  would  show  up  better.  In  a  code  in  which  there  were 
many  detectors  and  every  collision  was  used  to  obtain  an  estimate  at  every 
detector,  the  reverse  would  be  true.  We  believe  that  this  is  not  a  sensible 
thing  to  do;  the  next  section  discusses  the  alternatives. 

Unfortunately  it  is  not  possibie  to  use  a  simple  criterion  like  Q  to  compare 
estimation  method  1  with  the  others.  Nevertheless,  we  note  that  almost  all  the 
calculations  were  arranged  so  that  the  same  time  was  used  in  a  method  1  run  as 
in  a  method  2.  With  one  or  two  exceptions,  the  deviation  from  the  right  answer  is 
much  larger  in  1 .  As  we  have  seen,  the  infinite  variance  method  converges  more 
slowiy  so  that  additional  sampling  will  enlarge  the  difference. 

Fig.  3.1  shows  a  comparison  of  cumulative  averages  of  flux  at  a  point  as 
given  by  methods  1  and  2a.  The  results  are  for  a  problem  with  C  =  0.3  and  a  point 
0.5  mean  free  paths  from  the  source.  We  see  that  the  method  1  averages  are  con¬ 
sistently  poorer  by  far  and  furthermore  are  scarcely  better  after  many  thousand 
histories  than  after  a  few  hundred. 

On  this  basis  it  is  recommended  that  method  2a  or  2b  be  used  if  flux  at  a 
point  must  be  determined.  The  experiments  show  2a  better  than  2b,  at  least  in 
the  version  of  the  latter  that  was  used.  Until  a  better  method  along  the  lines  of 
2b  is  worked  out,  we  suggest  use  of  2a. 
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Fig.  3,1  —  Cumulative  average  collided  flux  after  n  histories 


4.  FLUX  ESTIMATION  AT  RANDOM 

In  calculating  the  radiation  field  of  a  weapon,  it  is  likely  that  answers  at 
many  detectors  would  be  required.  Suppose  that  a  flux  estimation  method  is  used 
that  gives  an  answer  for  every  collision.  It  is  then  possible,  or  even  likely,  that  a 
large  part  of  the  computation  time  will  be  taken  up  with  this  flux  estimation.  Much 
of  this  time  will  be  wasted,  since  any  given  collision  will  generally  be  far  from 
most  of  the  detectors.  The  partial  answers  at  these  detectors  will  be  very  small, 
and  the  computation  time  used  would  be  better  used  following  more  histories. 

It  appears  to  be  desirable,  then,  to  cut  off  the  flux  estimation  process  at 
large  distances.  One  way  of  doing  this  would  be  simply  not  to  carry  out  the  esti¬ 
mation  when  the  separation  of  collision  and  detector  is  larger  than  a  certain  num¬ 
ber  of  mean  free  paths.  This  results  in  a  consistent  underestimation  of  the  flux, 
but  if  the  cutoff  is  at  a  sufficiently  large  distance  the  error  becomes  small. 

Alternatively,  we  may  choose  to  carry  out  the  flux  estimation  with  a  proba¬ 
bility  which  is  a  function  of  collision-detector  separation.  That  is,  let  R  be  this 
separation  and  let  q(juR)  «  1  be  the  probabilityof  making  an  estimate.  If  q  becomes 
small  as  jiR  becomes  large,  then  a  cutoff  is  accomplished.  It  seems  plausible  that 
q  should  be  very  close  to  1  for  small  argument  (say  within  a  mean  free  path)  and 
decrease  monotonically. 

Suppose  that  when  a  particle  goes  into  collision  at  P  an  estimate  f(P)  is 
made  for  a  particular  flux.  Then  the  average  answer  is 

<f>  =  /  ^  (P)  f(P)  dP  (4.1) 
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If  we  carry  out  the  estimate  with  probability  of  q(P)  then  the  estimate  is  f(P)/q(P) 
so  that 

<f>  -  /  tjj  (P)  q(P)  [f(P)/q(P)]  dP.  (4.2) 

The  mean  square  answer  is 

<f2>  =  J  ip  (P)  -jH  dP  +  /  dP  ^(P)  f(P)  /  dP'  ^(P'  iP)  f(P0  (4.3) 

~  q(P)  ~  ~~~ 

where  tpiP'  IP)  is  the  collision  density  at  P'  given  that  a  collision  took  place  at  P. 
Clearly 

i^(P'IP)  =  C(E-E',n-n'lr)  T(r-rME',n')  (4.4) 

Eq.  4.3  shows  that  <f^>  and  hence  the  variance  in  one  history  are  smallest  if  q  =  1 
throughout.  On  the  other  hand  the  computing  time  is 

T  =  To+t,/  j/.(P)q(P)  dP  (4.5) 

where  Tq  is  the  time  for  tracing  the  histories  and  computing  and  testing  q.  Since 
the  time  decreases  as  q  is  made  smaller,  and  since  our  general  aim  is  to  mini¬ 
mize  variance  times  time,  it  is  possible  that  q  5^  1  is  optimum. 

The  determination  of  the  best  q  from  Eqs.  4.4  and  4.5  appears  very  difficult. 
In  the  brief  work  devoted  to  this  aspect  of  scoring,  we  have  used  simple  functional 
forms  for  q,  each  containing  an  adjustable  parameter. 

The  work  on  comparing  various  scoring  methods  for  flux  at  a  point  showed 
that  it  requires  very  lengthy  sampling  (at  least  on  a  Datatron)  to  compute  relative 
variances  accurately.  We  tried,  therefore,  to  study  the  effects  of  using  q  analyti¬ 
cally. 

For  a  straight  ahead  transport  problem  with  absorption,  the  flux  satisfies 


^  +  /id)  =  C  /ud) 
dx 


with  the  solution 


ip= 


-(l-C)/iX 
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We  may  use  ip  together  with  the  estimates  for  the  current  at  X 


and 

^-/3fi(X-x)  (4.6) 

in  order  that  the  integrals  in  Eqs.  4.3  and  4.5  may  be  done  analytically. 

In  this  way  it  is  possible  to  study  the  behavior  of  Q  (variance  x  time)  as  a 
function  of  C,  fjx,  and  /3. 

For  a  hypothetical  problem  with  many  detectors,  a  weighted  average  of  the 
individual  variances  is  used  instead  of  a  single  variance.  If  is  the  variance  at 
the  ith  detector  (at  X^)  then  the  efficiency  is  computed  from 


and  the  time  includes  flux  estimation  for  each  detector. 

The  results  will  not  be  shown  in  detail  since  they  are  not  extensive  enough 
for  interpretation  and  apply  to  a  model  probably  too  artificial  to  believe.  They 
indicate,  however,  that  savings  on  the  order  of  a  factor  of  1.4  to  5  might  be 
achieved.  The  larger  factor  is  associated  with  many  detectors. 

A  similar  calculation  was  carried  out  with  a  collision  density  altered  by 
importance  sampling.  The  general  conclusion  is  the  same:  by  proper  choice  of 
13,  savings  of  a  factor  of  up  to  about  10  might  be  expected. 

It  was  apparent  in  both  of  these  that  the  detailed  results  are  sensitive  to  the 
details  of  the  assumptions  of  distribution  of  collisions,  as  well  as  q  and  so  on. 
Apart  from  the  suggestion  that  the  method  is  worth  pursuing  further,  we  believe 
that  more  useful  conclusions  may  not  be  drawn. 

In  an  effort  to  apply  the  idea  to  a  more  realistic  case,  one  machine  code  for 
estimating  flux  at  a  point  was  modified  to  carry  out  the  scoring  at  random.  For 
this  work  the  function 
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(4.7) 


1  +  (niiRy 

was  used. 

This  q  has  two  advantages.  For  one,  it  falls  below  one  slowly  near  R  =  0. 
Secondly,  it  makes  the  computation  of  variance  as  a  function  of  easier.  From 
Eq.  4.3  we  see  that  the  variance  must  have  the  functional  form 

V  =  A  +  B/3^ 

By  plotting  experimental  variances  as  a  function  of  ^  it  is  possible  to  determine 
a  reasonable  straight  line  in  spite  of  fluctuations.  See  Fig.  4.1.  It  is  straightfor¬ 
ward  to  determine  the  computing  time  as  a  function  of  i3  as  in  Fig.  4.2.  The  times 
shown  refer  to  the  computation  for  a  single  detector,  but  the  average  time  for  col¬ 
lision  mechanics  may  be  found.  From  this  the  time  for  a  many  detector  problem 
(with  each  detector  the  same  distance  from  the  source)  may  be  computed.  The  Q 
is  shown  in  Figs.  4.3  and  4.4  for  one  and  10  detectors,  respectively,  all  two  mean 
free  paths  from  the  source.  The  savings  are  at  most  about  24%  which  is  not  very 
encouraging.  However,  since  the  penetrations  are  not  very  great,  not  much  should 
be  expected. 

The  potential  of  the  notion  of  doing  the  flux  calculation  at  random  remains  to 
be  proved.  It  would  be  best  to  try  the  experiments  in  a  code  nearly  like  the  final 
transport  code.  In  view  of  the  many  possibilities  of  choice  of  q,  it  would  be  very 
desirable  to  get  some  insight,  if  only  qualitative,  of  the  characteristics  of  an  opti¬ 
mum  q.  It  appears  that  the  method  is  likely  to  help  shorten  the  computation.* 


*Use  of  a  q  may  be  even  more  justified  in  problems  in  which  the  flux  estimation 
time  increases  with  collision-detector  separation.  Any  problem  with  complicated 
geometry  will  have  this  characteristic. 
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Relative  Variance 


0.01  0.05  0.1  0.15  0.2  0.25  0.3  0.35  0.4 


Fig.  4.1  —  Flux  estimation  at  random;  relative  variance  for  one 
history  vs 
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Q/<</7> 


/3 


Fig.  4.3 —  Flux  estimation  at  random;  variance  x  time  for  one 
detector 
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Fig.  4.4 —  Flux  estimation  at  random;  variance  x  time  for  10  de¬ 
tectors 
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